The present invention relates to the field of seismic data processing. In particular, this invention relates to a method of efficiently and accurately producing the filter coefficients for wavefield extrapolation particularly useful for migrating seismic data using a multiprocessor parallel computer.
1. Seismic Acquisition
The Earth's subsurface can be imaged by a seismic survey, therefore, seismic data acquisition and processing are key components in geophysical exploration. In a seismic survey, elastic acoustic waves are generated by a source at the Earth's surface and the waves are radiated into the Earth's subsurface. For land seismic surveys, the usual source is dynamite or a seismic vibrator, while for a marine seismic survey the source is typically an airgun array.
As the waves radiate downward through the Earth's subsurface, they reflect and propagate upwards towards the surface whenever the subsurface medium changes. The upward reflections are detected by a number of receivers and the reflected data recorded and processed in order to image the subsurface. Interpretation of these acoustic images of the subsurface formation leads to the structural description of the subsurface geological features, such as faults, salt domes, anticlines, or other features indicative of hydrocarbon traps.
While two dimensional ("2D") seismic surveys have been conducted since the 1920's, three dimensional ("3D") seismic surveys have only recently become widely used. 3D surveys more accurately reflect the subsurface positions of the hydrocarbon traps, but are expensive and time consuming to acquire and process. For an offshore 3D data set covering a 20.times.20 km area, it costs about $3M dollars ( 1991 dollars) to acquire the data with another $1 M dollars for data processing to transform the raw data into usable images. Because the cost of such a seismic survey is considerably less than the cost of drilling an offshore oil well, 3D seismic surveys are often worth the investment.
Although 3D marine surveys vary widely in size (1,000 to 100,000 km.sup.2), a typical marine survey might generate in excess of 40,000 data acquisition tapes. Data is accumulated at a staggering rate, about 1.5 million data samples every 10 seconds. A significant amount of time and money is spent in processing this enormous amount of data.
The result of the seismic survey is thus an enormous amount of raw data indicative of reflected signals which are a function of travel time, propagation, and reflection affects. The goal is to present the reflected amplitudes as a function of lateral position and depth.
2. Seismic Processing
A typical marine seismic survey goes through three distinct sequential stages--data acquisition, data processing, and data interpretation. Data processing is by far the most time consuming process of the three. The acquisition time for a medium to large 3D marine seismic survey is in the order of two months. Data is acquired by survey vessels traversing an area of the ocean along a series of parallel lines. A vessel may tow a number of sources (usually airgun arrays) and a number of receiver strings called hydrophone streamers (of length up to 5 kilometers). Sources are fired at 5 to 10 second intervals and the reflected seismic waves measured by up to 1000 hydrophone groups in every streamer. The measurements are recorded digitally on magnetic tapes. In addition to seismic data, navigation information is also recorded for accurate positioning of the sources and receivers. The resulting digital data must then be rendered suitable for interpretation purposes by processing the data at an onshore processing center. The processing sequence can be divided into the following five processing steps.
1. Quality Control, filtering and deconvolution. This processing is applied on a trace basis to filter noise. sharpen the recorded response, suppress multiple echoes, and generally improve the signal-to-noise ratio. Most of these signal processing operations can be highly vectorized. PA0 2. Velocity analyses for migration. This processing estimates the velocity of the subsurface formations from the recorded data by modeling the propagation of acoustic waves with estimated velocities and checking for signal coherence in the acquired data. It is similar to migration but is applied to a small section of the data cube. PA0 3. D dip moveout correction and stacking. This processing step, generally the most input/output intensive part of the processing, (i) sums together several traces in order to eliminate redundancy and reduce the signal-to-noise ratio, (ii) corrects for time delays that occur when the reflected signal is recorded by successive hydrophones that are located increasingly farther away from the energy source, and (iii) positions and orients the stacked data in accordance with the navigation information. After this processing step, the data is referred to as stacked data. This step normally constitutes on the order of a 100 to 1 reduction in data volume. PA0 4. Migration. This processing step, computationally the most intensive, relocates the position of reflected strata, that are recorded in time, to their correct position in depth. PA0 5. Enhancement and filtering. This processing step is used to enhance the migrated data using digital Flitering techniques.
The stacking process (step 3) reduces the amount of data to what is essentially a three dimensional array of numbers (i.e. a data cube) representing amplitudes of reflected seismic waves recorded over a period of time (usually 8 seconds). Such data cubes can be large, for example, a medium size 3D survey may produce cubes as large as 1000.times.1000.times.2000 of floating-point numbers.
The stacked data cube represents a surface recording of acoustic echoes returned from the earth interior and is not usually directly interpretable. The migration (or acoustic imaging process, step 4) is used to convert stacked data into an image or a map which can then be viewed as a true depth map cut out of the survey area.
Thus, migration is one of the most critical and most time consuming components in seismic processing is migration. Generally speaking, migration transforms the seismic data recorded as a function of time into data positioned as a function of depth using preliminary knowledge of the propagation velocities of the subsurface. In particular, migration moves dipping reflectors to their true subsurface position. Migration is typically performed on post stack seismic data to reduce the amount of processing time, but even so takes weeks of conventional supercomputer time for even medium size post stack seismic data cubes.
Most of the migration methods are based on the one way acoustic wave equation (compressional waves considered, shear waves ignored) using the exploding reflector model. In the exploding reflector model, stacked data are assumed to be recordings of a multitude of sources distributed along geological boundaries and exploded simultaneously. The post stack seismic data cube is considered to be recordings of upward traveling waves as they emerge from the Earth. (See generally, J. F. Claerbout, Imaging the Earth's Interior (1985), M. Dobrin & C. Savit, Geophysical Prospecting (1988), R. E. Sheriff, Geophysical Methods (1989), incorporated by reference for background).
3. Wavefield Extrapolation
Wavefield extrapolation is critical in the migration process. (The terms "wavefield extrapolation" and "downward continuation" are sometimes used interchangeably.) In wavefield extrapolation an imaging step is performed at each depth, which is the extraction of zero time amplitudes from the downward continued data cube. Wavefield extrapolation in the space-frequency (x,y,f) domain is preferable because it can more accurately image steeply dipping geological layers and certain mathematical operations are easier. In the space-frequency domain, the 2D scalar (two way) wave equation can be written as: ##EQU1##
wherein k.sub.2 =(w.sup.2 /c.sup.2 -k.sub.x.sup.2).sup.1/2, U=U (w,k.sub.x,z), D=(w,k.sub.x,z) represent the upgoing and downgoing waves respectively, w is the frequency (measured in radians per unit time), c is propagation velocity, and k.sub.x is the wavenumber (measured in radians per sample) in the x direction. r(k.sub.z,z) is the reflectivity function. This equation holds for a horizontal layered media. The first term in the equation accounts for the one way propagation of a wave in a homogeneous media. The second term accounts for transmission losses and coupling between the upgoing and downgoing waves at the interfaces. If transmission losses are neglected, the 2D scalar wave equation can be rewritten as: EQU .differential.P/.differential.z=.+-.ik.sub.z P
where P may be upward U or downward D. This is the basis for one way wave propagation such as using the exploding reflector model.
The analytical solution to this equation is: EQU P(w,k.sub.x,z+.DELTA.z)=e.sup..+-.ik.sup..sub.z .sup..DELTA.z P(w,k.sub.x,z)
corresponding to downward extrapolation of one way waves. This analytical solution can be approximated in the space-frequency domain with known finite difference techniques. These finite difference techniques resemble a 2D convolution when known techniques have been applied, collectively known as splitting or layer splitting. However, using conventional layer splitting techniques, computer accuracy or computer efficiency is sacrificed. It is useful to use a finite difference approximation in the space-frequency domain and recasts the problem as a filter. The Fourier transform approximates: ##EQU2##
where w denotes temporal frequency, .nu. is the velocity, and z and x are vertical and horizontal spatial sampling intervals, and k is the wavenumber frequency. Wavefield extrapolation may be reduced to applying a filter with the characteristics of equation (1) independently to every frequency in the data cube.
Implicit filtering methods are widely used to extrapolate seismic wavefields in depth because of their stability. Such implicit methods for depth extrapolation do not permit the amplitude of the extrapolated wavefield to grow with depth. However, explicit filtering methods can be more easily implemented and are believed to be more efficiently executed on different computer architectures. In addition to simplicity, explicit methods are more readily extendable to 3D depth migration.
Because of the computational advantages of explicit methods, different methods have been proposed. O. Holberg, Towards Optimum One-way Wave Propagation, 306 Geophysical Prospectizg 99-114 (1988) generally discusses the wavefield extrapolation problem for depth migration, pre or post stack. Holberg reviews the various finite difference explicit methods for approximating the acoustic one-way equation in the space-frequency domain. Holberg then proposes a least-squares solution which minimizes the total phase and amplitude errors and their gradients for a range of propagation angles.
Building on Holberg's work, D. Hale Stable Explicit Depth Extrapolation of Seismic Wavefields, 56 Geophysics 1770 (1991), proposed a modified Taylor series expansion method for defining the coefficients to approximate the extrapolation filter of equation(1). Hale's modified Taylor series solution apparently yields a stable extrapolation filter approximation. Hale claims several advantages of the modified Taylor series method over the Holberg least-squares method. While the modified Taylor series method is stable, the phase error is increased. To reduce the phase error, Hale proposes increasing the number of coefficients in the extrapolation filter.
4. 3D Seismic Migration
Pieprzak and Highnam in U.S. patent application Ser. Nos. 07/811,414 (pending) and 07/811,565 (allowed) (incorporated by reference) describe accurate 3D one pass migration methods which employ the parallel features of a multiprocessor parallel computer to perform seismic migration. The 3D migration method employed by Pieprzak and Highnam. building on the work of G. Blacquiere, D. Hale, and J. McClellan, uses only the one dimensional extrapolation filters, so that the computational cost of 3D depth migration increases only linearly with N, the number of filter coefficients. The Pieprzak and Highnam migration method implements the use of a recursive Chebyshev filter structure for 2D convolution with radically symmetric operators using the 1D extrapolation filters. The method employed by Pieprzak and Highnam, both for time and depth, employs conventional wavefield extrapolation techniques, such as Holberg's least squares method, to compute the filter coefficients used in its two dimensional recursive filter approximation. The Pieprzak and Highnam migration runtime is totally dominated by the time necessary to apply these filter coefficients. Because of such extensive use, improvements in the accuracy of these filter coefficients increases the accuracy of the migration.
It is, therefore, an object of the present invention to provide an efficient method of approximating such wavefield extrapolation filter coefficients to increase computational efficiency, accuracy, or both. It is another object of the invention to automatically provide extrapolation filter coefficients for use in seismic migration. Still another object of the invention is to provide extrapolation filters useful for 3D seismic migration using a multiprocessor parallel computer.
In this application all references to patent documents are incorporated by reference and all references to non patent documents are incorporated by reference for background.